Análisis Exploratorio de Datos

Cargar paquetes

library(readxl)
library(tidyverse)
library(kableExtra)
library(sf)
library(mapview)
library(vegan)
library(devtools)
library(indicspecies)
library(SpadeR)
library(adespatial) #TBI, para comparar matrices de comunidad entre estaciones
library(RColorBrewer)
library(spdep)
download.file('https://raw.githubusercontent.com/biogeografia-master/scripts-de-analisis-BCI/master/biodata/funciones.R', 'funciones1.R', )
source('funciones1.R')
source('funciones2.R')
source('estadistica-zonal.R')

Cargar datos

tabla_odk <- read_csv('tabla ODK.csv')
tabla_odk <- tabla_odk[grep('prueba', tabla_odk$codigo, ignore.case = T, invert = T), ]
tabla_odk$codigo[grep('h14m75', tabla_odk$codigo, ignore.case = T)] <- 'H10M75'
tabla_odk$codigo[grep('h1m14', tabla_odk$codigo, ignore.case = T)] <- 'H12M14'
tabla_odk$codigo <- toupper(tabla_odk$codigo)
tabla_odk$codigo[grep('h1m88', tabla_odk$codigo, ignore.case = T)] <- 'H2M88'
tabla_odk$codigo[grep('2h23h42', tabla_odk$codigo, ignore.case = T)] <- '2H23M42'
tabla_odk$codigo[grep('2h23h43', tabla_odk$codigo, ignore.case = T)] <- '2H23M43'
tabla_odk$codigo[grep('2H25M147', tabla_odk$codigo, ignore.case = T)] <- '2H24M147'

datos_orig <- read_xlsx(
  path = 'Datos de muestreo Marzo y Mayo 2022.xlsx',
  col_types = c('numeric', rep('text', 7), 'date', 'numeric', rep('text', 5)))
which(is.na(match(tabla_odk$codigo, datos_orig$Código))) #Todos los elementos tienen contrapartes entre tabla si "integer(0)"

## integer(0)

ncol(datos_orig) # Número de columnas

## [1] 15

nrow(datos_orig)

## [1] 255

colnames(datos_orig)

##  [1] "Parcela"          "Código"           "Coordenadas"      "Datum"           
##  [5] "Altura"           "Incertidumbre"    "Especie"          "Colector/os"     
##  [9] "Fecha"            "No. de individos" "Observaciones"    "País"            
## [13] "Provincia"        "Municipio"        "Área protegida"

datos_01 <- datos_orig %>%
  mutate(lat = str_split(datos_orig$Coordenadas, pattern = ' ', simplify = T)[,1],
         lon = str_split(datos_orig$Coordenadas, pattern = ' ', simplify = T)[,2]) %>% 
  mutate_at(c('lat', 'lon'), as.numeric)
datos_02 <- datos_01 %>%
  select(-Parcela) %>%
  mutate(Hexágono = as.numeric(gsub('.*H', '', gsub('M.*', '', Código)))) %>% 
  select(Hexágono, everything())
datos_02 %>%
  select(Hexágono, Código, Especie) %>% 
  head() %>% 
  kable()

Hexágono	Código	Especie
13	H13M1	Pheidole subarmata
13	H13M2	Dorymyrmex antillanus
13	H13M3	Paratrechina longicornis
13	H13M4	Pheidole subarmata
13	H13M5	Solenopsis geminata
13	H13M6	Pheidole jelskii

datos_03 <- datos_02 %>%
  inner_join(
    tabla_odk %>%
      select(codigo, `coordenadaid-Latitude`, `coordenadaid-Longitude`), by = c('Código' = 'codigo')) %>% 
  select(-lat,-lon) %>% 
  rename(lat = `coordenadaid-Latitude`,
         lon = `coordenadaid-Longitude`)
datos_03[grep('^H13M1$', datos_03$Código, ignore.case = T), 'Hexágono'] <- 12 #Misplaced, corrected
datos_03$Código[grep('^H13M1$', datos_03$Código, ignore.case = T)] <- 'H12M1' #Misplaced, corrected
datos_03[grep('^H1M87$', datos_03$Código, ignore.case = T), 'Hexágono'] <- 4 #Misplaced, corrected
datos_03$Código[grep('^H1M87$', datos_03$Código, ignore.case = T)] <- 'H4M87' #Misplaced, corrected
datos_03$Especie <- gsub('Pheidole  jelskii', 'Pheidole jelskii', datos_03$Especie)
datos_03$Especie <- gsub('Ontomachus ruginosa', 'Odontomachus ruginodis', datos_03$Especie)
datos_03$Especie <- gsub('Monomorium floricula', 'Monomorium floricola', datos_03$Especie)
datos_03$Especie <- gsub('Monomorium lanuginosa', 'Tetramorium lanuginosum', datos_03$Especie)
# datos_03 %>% write_csv('datos_03.csv')
# write_excel_csv(datos_03, 'datos_03_compatible_con_excel.csv')

Crear mapa de las colectas

# datos_03 <- read_csv('datos_03.csv')
datos_sf <- datos_03 %>% 
  st_as_sf(coords = c('lon', 'lat'))
st_crs(datos_sf) <- 4326
datos_sf %>% plot()

# st_write(datos_sf, 'datos_sf.gpkg', delete_dsn = T)

# Generar mapa de cuadros sin simbología
## Centroid
# datos_sf %>% mutate(i=1) %>% 
#     group_by(i) %>% 
#     summarize(geometry = st_union(geometry)) %>% 
#     st_centroid %>% mutate(lng )
## Mapa
mapa_cuadros <- mapView(
  datos_sf,
  col.regions = 'grey80',
  alpha.regions = 0.3,
  map.types = 'Esri.WorldImagery',
  legend = F, zoom = 14,
  zcol = 'Código') %>% leafem::addStaticLabels(label = datos_sf$Código) # %>%
  # leaflet::setView(
  #   lng = -79.85136,
  #   lat = 9.15097,
  #   zoom = 15)
mapa_cuadros

# mapa_cuadros %>% mapshot(file = 'mapa_punto.png') #Genera archivo
# Capa de hexagonos
hex <- st_read('cuadricula-final.gpkg')

## Reading layer `cuadricula-final' from data source 
##   `C:\Users\User\Desktop\Universidad\Tesis\Datos de los muestreos\cuadricula-final.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 28 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 402697.2 ymin: 2041872 xmax: 403143.5 ymax: 2042260
## Projected CRS: WGS 84 / UTM zone 19N

hex %>% inner_join(#NUMERO DE NIDOS POR HEXAGONO
  datos_03 %>%
    group_by(Hexágono) %>% 
  summarise(n = n()), by = c('ENLACE' = 'Hexágono')) %>% 
  mapview(zcol = 'n')

Tabla de número de nidos por Hexágono

datos_03 %>% group_by(Hexágono) %>%
  summarise(n = n()) %>% arrange(desc(n)) %>% 
  head() %>% 
  kable()

Hexágono	n
12	19
16	16
21	16
20	15
13	13
6	12

Crear matriz de comunidad

# datos_03 <- read_csv('datos_03.csv') #Descomentar si se desea comenzar el script desde este punto (no olvidar cargar paquetes)
datos_03

## # A tibble: 255 × 17
##    Hexágono Código Coordenadas            Datum Altura  Incert…¹ Especie Colec…²
##       <dbl> <chr>  <chr>                  <chr> <chr>   <chr>    <chr>   <chr>  
##  1       12 H12M1  18.4675906 -69.9139862 WGS84 53 msnm 5 metros Pheido… D. Guz…
##  2       13 H13M2  18.4672803 -69.9192371 WGS84 54 msnm 4 metros Dorymy… D. Guz…
##  3       13 H13M3  18.4672054 -69.9192181 WGS84 54 msnm 5 metros Paratr… D. Guz…
##  4       13 H13M4  18.4669508 -69.9193136 WGS84 54 msnm 7 metros Pheido… D. Guz…
##  5       13 H13M5  18.4671236 -69.9191236 WGS84 54 msnm 3 metros Soleno… D. Guz…
##  6       13 H13M6  18.4671876 -69.9193009 WGS84 53 msnm 4 metros Pheido… D. Guz…
##  7       12 H12M7  18.4675406 -69.9194121 WGS84 53 msnm 7 metros Pheido… D. Guz…
##  8       12 H12M8  18.4675806 -69.9193971 WGS84 53 msnm 5 metros Pheido… D. Guz…
##  9       12 H12M9  18.4675541 -69.9193724 WGS84 53 msnm 7 metros Pheido… D. Guz…
## 10       12 H12M10 18.4678106 -69.9148187 WGS84 53 msnm 7 metros Soleno… D. Guz…
## # … with 245 more rows, 9 more variables: Fecha <dttm>,
## #   `No. de individos` <dbl>, Observaciones <chr>, País <chr>, Provincia <chr>,
## #   Municipio <chr>, `Área protegida` <chr>, lat <dbl>, lon <dbl>, and
## #   abbreviated variable names ¹​Incertidumbre, ²​`Colector/os`

generar_mc <- function(patron = '^H.*', presencia_ausencia = T){
  datos_03 %>%
    filter(grepl(patron, Código)) %>% 
    select(Hexágono, Especie) %>%
    mutate(n = 1) %>% 
    pivot_wider(names_from = Especie, values_from = n,
                values_fn = ~ sum(.x, na.rm = TRUE),
                values_fill = 0) %>%
    arrange(Hexágono) %>% 
    column_to_rownames('Hexágono') %>% 
    {if(presencia_ausencia) replace(., . > 1, 1) else .}
}
mc_marzo <- generar_mc(patron = '^H.*')
# png('grafico_de_mosaicos_marzo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(
  add_row(mc_marzo, .before = 1) %>% replace(is.na(.), 0) %>% remove_rownames()
)

# dev.off()
mc_mayo <- generar_mc(patron = '^2H.*')
# png('grafico_de_mosaicos_mayo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_mayo)

# dev.off()
mc_pooled <- generar_mc(patron = '.*')
# png('grafico_de_mosaicos_pooled.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_pooled)

# dev.off()
mc_marzo_nidos <- generar_mc(patron = '^H.*', presencia_ausencia = F)
# png('grafico_de_mosaicos_nidos_marzo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(
  add_row(mc_marzo_nidos, .before = 1) %>% replace(is.na(.), 0) %>% remove_rownames()
)

# dev.off()

mc_mayo_nidos <- generar_mc(patron = '^2H.*', presencia_ausencia = F)
# png('grafico_de_mosaicos_nidos_mayo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_mayo_nidos)

# dev.off()

# Nidos combinados
mc_nidos <- generar_mc(patron = '.*', presencia_ausencia = F)
# png('grafico_de_mosaicos_nidos.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_nidos)

# dev.off()

# Reescalado
mc_puente <- scales::rescale(as.matrix(mc_nidos), to= c(1, 5))
mc_nidos_reescalado <- round(ifelse(mc_puente==1, 0, mc_puente), 0)

EDA

mc_todas <- list(mc_marzo, mc_mayo, mc_pooled)
names(mc_todas) <- c('Marzo', 'Mayo', 'Combinada')
#' ### Lista de especies
map(mc_todas, 
    ~ sort(colnames(.x)) %>%
      as_tibble() %>%
      mutate(fila = 1:nrow(.)) %>%
      select(2, 1) %>% 
      kable(col.names = c('#', 'Especie'), table.attr='class="myTable"'))

$Marzo

#	Especie
1	Cyphomyrmex rimosus
2	Dorymyrmex antillanus
3	Mycetomoellerius jamaicensis
4	Odontomachus ruginodis
5	Paratrechina longicornis
6	Pheidole jelskii
7	Pheidole moerens
8	Pheidole subarmata
9	Solenopsis geminata
10	Tetramorium bicarinatum
11	Tetramorium lanuginosum
12	Wasmannia auropunctata

$Mayo

#	Especie
1	Acropyga dubitata
2	Brachymyrmex heeri
3	Cyphomyrmex rimosus
4	Dorymyrmex antillanus
5	Monomorium floricola
6	Mycetomoellerius jamaicensis
7	Odontomachus ruginodis
8	Paratrechina longicornis
9	Pheidole jelskii
10	Pheidole subarmata
11	Pogonomyrmex schmitti
12	Solenopsis geminata
13	Tetramorium lanuginosum
14	Wasmannia auropunctata

$Combinada

#	Especie
1	Acropyga dubitata
2	Brachymyrmex heeri
3	Cyphomyrmex rimosus
4	Dorymyrmex antillanus
5	Monomorium floricola
6	Mycetomoellerius jamaicensis
7	Odontomachus ruginodis
8	Paratrechina longicornis
9	Pheidole jelskii
10	Pheidole moerens
11	Pheidole subarmata
12	Pogonomyrmex schmitti
13	Solenopsis geminata
14	Tetramorium bicarinatum
15	Tetramorium lanuginosum
16	Wasmannia auropunctata

Matriz ambiental

cobertura <- st_read('coberturas.gpkg')

## Reading layer `coberturas' from data source 
##   `C:\Users\User\Desktop\Universidad\Tesis\Datos de los muestreos\coberturas.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 82 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 402699.2 ymin: 2041878 xmax: 403116 ymax: 2042246
## Projected CRS: WGS 84 / UTM zone 19N

hex_inters <- st_as_sf(st_intersection(st_union(cobertura), hex)) %>% mutate(ENLACE = hex$ENLACE)
ma_sf <- zstatsf(zones = hex_inters, values = cobertura, grpx = "ENLACE", grpy = "tipo")
ma_sf <- ma_sf %>% mutate(across(2:5, units::drop_units))
ma_sf %>% plot()

ma <- ma_sf %>% st_drop_geometry() %>% column_to_rownames('ENLACE')
ma

##    construido, mobiliario (bordes edificios, acerado, bancos, postes...)
## 1                                                              0.0000000
## 2                                                             18.4869756
## 3                                                              0.0000000
## 4                                                              3.4888882
## 5                                                              0.0000000
## 6                                                              2.4769313
## 7                                                              0.6085964
## 8                                                             27.2501334
## 9                                                              7.5820449
## 10                                                             0.1769173
## 11                                                             5.1515836
## 12                                                            38.0738217
## 13                                                            51.8993681
## 14                                                             0.0000000
## 15                                                             1.1269871
## 16                                                             3.7630314
## 17                                                            69.1243002
## 18                                                            16.2492695
## 19                                                            25.6929698
## 20                                                            15.8724897
## 21                                                            21.0199767
## 22                                                            43.0273268
## 23                                                            11.3836316
## 24                                                            44.3322007
## 25                                                            23.9599716
## 26                                                             1.7747530
## 27                                                            18.4280883
## 28                                                             0.3534807
##       dosel edificación erguida suelo, herbáceas, no edificado ni cubierto
## 1  89.86997            0.000000                                10.13003022
## 2  80.17277            0.000000                                 1.34025095
## 3  99.88796            0.000000                                 0.11204383
## 4  61.57379           24.838187                                10.09913880
## 5  77.72845           22.271546                                 0.00000000
## 6  95.53464            0.000000                                 1.98843129
## 7  97.26926            0.000000                                 2.12214164
## 8  67.74032            0.000000                                 5.00955070
## 9  54.37792           37.902586                                 0.13744975
## 10 88.68445            0.000000                                11.13863183
## 11 80.03134            1.765534                                13.05153872
## 12 52.72902            0.000000                                 9.19716130
## 13 46.72021            1.068254                                 0.31216408
## 14 96.50561            0.000000                                 3.49438827
## 15 97.93200            0.000000                                 0.94101640
## 16 96.23697            0.000000                                 0.00000000
## 17 26.98363            0.000000                                 3.89206487
## 18 63.99899           19.294636                                 0.45710407
## 19 74.29148            0.000000                                 0.01555076
## 20 84.12751            0.000000                                 0.00000000
## 21 67.26251            0.000000                                11.71751267
## 22 15.30474           28.156893                                13.51103755
## 23 41.44998           46.221870                                 0.94452067
## 24 53.38524            0.000000                                 2.28256075
## 25 76.04003            0.000000                                 0.00000000
## 26 65.53121           24.489490                                 8.20454626
## 27 78.54138            0.000000                                 3.03053153
## 28 93.44066            0.000000                                 6.20586141

rowSums(ma) #Comprobación de suma 100 por filas

##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
## 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 
##  21  22  23  24  25  26  27  28 
## 100 100 100 100 100 100 100 100

Matriz combinada

Resúmenes

#' ### Número de sitios, tanto en matriz de comunidad como en ambiental
#' Verifica que coinciden
nrow(mc_todas$Combinada) #En la matriz de comunidad

## [1] 28

nrow(ma) #En la matriz ambiental

## [1] 28

#' ### Riqueza numérica de especies (usando matriz de comunidad) por quadrat
#' Nota: cargar paquete vegan arriba, en el área de paquetes
specnumber(mc_todas$Combinada)

##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
##  2  3  2  4  4  7  9  3  4  2  6  5  6  2  6  3  1  3  4  4  6  3  2  4  5  2 
## 27 28 
##  4  4

sort(specnumber(mc_todas$Combinada)) # Ordenados ascendentemente

## 17  1  3 10 14 23 26  2  8 16 18 22  4  5  9 19 20 24 27 28 12 25 11 13 15 21 
##  1  2  2  2  2  2  2  3  3  3  3  3  4  4  4  4  4  4  4  4  5  5  6  6  6  6 
##  6  7 
##  7  9

summary(specnumber(mc_todas$Combinada)) # Resumen estadístico

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.750   4.000   3.929   5.000   9.000

#' ### Riqueza numérica de toda la "comunidad"
specnumber(colSums(mc_todas$Combinada))

## [1] 16

Agrupamiento, asociacion de especies con habitat (Objetivo 3 o 4)

En esta seccion se realiza una clasificacion no supervisada de habitats en funcion de su composicion de especies de hormigas, y se explora cuales atributos (de los 4 generados a partir del mapa de cobertura) caracterizan a dichos habitats.

Filtrado para excluir hexagonos muy pobres o especies muy raras, si aplica

mc_orig <- mc_todas$Combinada
data.frame(`Número de hexágonos` = sort(colSums(mc_orig), decreasing = T), check.names = F) %>% 
  kable(booktabs=T) %>%
  kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
  gsub(' NA ', '', .) # Número de hexágonos en los que está presente cada especie

	Número de hexágonos
Pheidole subarmata	25
Solenopsis geminata	20
Dorymyrmex antillanus	19
Paratrechina longicornis	8
Pheidole jelskii	8
Wasmannia auropunctata	6
Cyphomyrmex rimosus	6
Odontomachus ruginodis	5
Tetramorium lanuginosum	3
Pogonomyrmex schmitti	3
Mycetomoellerius jamaicensis	2
Pheidole moerens	1
Tetramorium bicarinatum	1
Brachymyrmex heeri	1
Monomorium floricola	1
Acropyga dubitata	1

# Usa el vector anterior para determinar un umbral o rango de registros para filtrar tu matriz
# ¿En cuántos hexágonos está cada especie? Filtra tus datos usando tu propio criterio.
# Especies que aparecen en pocos hexágonos se consideran "raras". Por ejemplo, si una especie sólo
# aparece en un hexágono en todo el país, es un "singleton", si en dos, "doubleton", y así.
# Estas especies podrían contribuir a generar "ruido" en análisis posteriores, se recomienda excluirlas.
# Elige un valor mínimo (representado por único número entero) o por un rango de enteros (e.g. de 10 a 20),
# para seleccionar las especies que estén mejor representadas de acuerdo a tu criterio.
# Por ejemplo, si usas el valor m, el script considerará a este valor como "el número mínimo de hexágonos
# en los que está representada una especie, y creará una matriz de comunidad de especies seleccionadas
# que están presentes en m hexágonos o más. Si eliges un rango, por ejemplo [m,n], el script generará
# una matriz de comunidad que representadas un mínimo de m hexágonos y un máximo de n hexágonos.
# (ambos extremos inclusive).
en_cuantos_hex <- 2 #!!!ESTE UMBRAL EXCLUYE A ESPECIES QUE APARECEN EN SOLO 1 HEXAGONO, Y CONSERVA SOLO AQUELLAS QUE APARECEN EN 2 O MAS HEXAGONOS
# Explicación: "en_cuantos_hex <- X", donde X es el número de hexágonos mínimo donde cada especie
# debe estar presente. IMPORTANTE: elige TU PROPIO umbral.
{if(length(en_cuantos_hex)==1) selector <- en_cuantos_hex:max(colSums(mc_orig)) else
  if(length(en_cuantos_hex)==2)
    selector <- min(en_cuantos_hex):max(en_cuantos_hex) else
      stop('Debes indicar uno o dos valores numéricos')}
selector

##  [1]  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

mc_orig_seleccionadas <- mc_orig[, colSums(mc_orig) %in% selector]

# Mínimo número de especies por hexágono
data.frame(`Número de especies por hexágono` = sort(rowSums(mc_orig), decreasing = T), check.names = F) %>% 
  kable(booktabs=T) %>%
  kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
  gsub(' NA ', '', .) # Número de hexágonos en los que está presente cada especie

	Número de especies por hexágono
7	9
6	7
11	6
13	6
15	6
21	6
12	5
25	5
4	4
5	4
9	4
19	4
20	4
24	4
27	4
28	4
2	3
8	3
16	3
18	3
22	3
1	2
3	2
10	2
14	2
23	2
26	2
17	1

min_especies_por_hex <- 2 #!!!ESTA LINEA SACA LOS HEXAGONOS QUE TIENEN SOLO 1 ESPECIE, Y CONSERVA SOLO AQUELLOS QUE TIENEN 2 O MAS ESPECIES PRESENTES (EN EL CASO DE LA MATRIZ DE ESPECIES COMBINADAS, EXCLUYE EL HEXAGONO 17)
# Explicación: "min_especies_por_hex <- Y", donde Y es el número mínimo (inclusive) de especies
# que debe existir en cada hexágono. Por debajo de dicho valor, el hexágono es excluido.
mi_fam <- mc_orig_seleccionadas[rowSums(mc_orig_seleccionadas)>=min_especies_por_hex, ]
nrow(mi_fam)

## [1] 27

# mi_fam <- mc_orig_seleccionadas[!rowSums(mc_orig_seleccionadas)==0, ] #Elimina filas sin registros
# rowSums(mi_fam) #Riqueza por hexágonos con especies seleccionadas. Comentado por extenso
all(rowSums(mi_fam)>0) #Debe ser TRUE: todos los hexágonos tienen al menos 1 registro

## [1] TRUE

ncol(mi_fam) #Riqueza de especies

## [1] 11

# IMPORTANTE, quitar Formicidae. ¡¡Muchos nombres que habría que resolver si se quisiera publicar!!
# Usar nombres cortos o abreviados para las especies
# nombres_largos <- colnames(mi_fam)
# (colnames(mi_fam) <- make.cepnames(word(colnames(mi_fam), 1, 2)))
# (df_equivalencias <- data.frame(
#   nombre_original = nombres_largos,
#   abreviado = colnames(mi_fam)))

Análisis de agrupamiento.

En esta subseccion se realiza propiamente el analisis de agrupamiento segun cuatro metodos distintos (single, completate, upgma, Ward); se elige, al final, el metodo Ward. Se contrastan los metodos con remuestreo multiescalas por bootstrap.

mi_fam_t <- decostand(mi_fam, 'hellinger') #Hellinger
env <- ma[match(rownames(mi_fam), rownames(ma)), ]
all(rownames(mi_fam) == rownames(env)) #Si es TRUE, sigue adelante

## [1] TRUE

library(broom)
library(cluster)
library(gclus)
library(pvclust)
mi_fam_d <- vegdist(mi_fam_t, "euc")
lista_cl <- list(
        cl_single = hclust(mi_fam_d, method = 'single'),
        cl_complete = hclust(mi_fam_d, method = 'complete'),
        cl_upgma = hclust(mi_fam_d, method = 'average'),
        cl_ward = hclust(mi_fam_d, method = 'ward.D2')
)
par(mfrow = c(2,2))
invisible(map(names(lista_cl), function(x) plot(lista_cl[[x]], main = x, hang = -1, cex = 1)))

par(mfrow = c(1,1))
cutree(lista_cl$cl_ward, 2)

##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 18 19 20 21 22 23 24 25 26 27 
##  1  1  2  1  1  2  2  2  2  1  2  1  1  1  2  1  1  1  1  1  1  1  1  2  2  1 
## 28 
##  2

Remuestreo multiescalar por bootstrap

if(interactive()) dev.new()
cl_pvclust_ward <-
        pvclust(t(mi_fam_t),
                method.hclust = "ward.D2",
                method.dist = "euc",
                iseed = 999, # Resultado reproducible
                parallel = TRUE)

## Creating a temporary cluster...done:
## socket cluster with 3 nodes on host 'localhost'
## Multiscale bootstrap... Done.

# Añadir los valores de p
plot(cl_pvclust_ward, hang = -1)
# Añadir rectángulos a los grupos significativos
lines(cl_pvclust_ward)
pvrect(cl_pvclust_ward, alpha = 0.91, border = 4)

Guardando agrupamientos a archivos

w_dend_reord <- reorder.hclust(lista_cl$cl_ward, mi_fam_d)
grupos_ward_k2 <- as.factor(cutree(lista_cl$cl_ward, k = 2))
table(grupos_ward_k2) #¿Cuántos hexágonos hay en cada grupo?

## grupos_ward_k2
##  1  2 
## 17 10

plot(w_dend_reord, hang = -1)
rect.hclust(tree = w_dend_reord, k = 2)

grupos_ward_k3 <- as.factor(cutree(lista_cl$cl_ward, k = 3))
table(grupos_ward_k3) #¿Cuántos hexágonos hay en cada grupo?

## grupos_ward_k3
##  1  2  3 
##  6 10 11

plot(w_dend_reord, hang = -1)
rect.hclust(tree = w_dend_reord, k = 3)

# Guardaré estos vectores en archivos para reutilizarlos en *scripts* posteriores: 
# saveRDS(grupos_ward_k2, 'grupos_ward_k2.RDS')
# saveRDS(grupos_ward_k3, 'grupos_ward_k3.RDS') #No hubo patrones significativos de asociacion

Caractersticas ambientales de grupos, asociacion con habitat

ma_para_combinadas <- ma[match(rownames(mi_fam), rownames(ma)), ]
(m_amb_ward_k2 <- ma_para_combinadas %>%
   rownames_to_column('hex_id') %>% 
   mutate(grupos_ward_k2) %>%
   pivot_longer(-c(grupos_ward_k2, hex_id), names_to = "variable", values_to = "valor"))

## # A tibble: 108 × 4
##    hex_id grupos_ward_k2 variable                                          valor
##    <chr>  <fct>          <chr>                                             <dbl>
##  1 1      1              construido, mobiliario (bordes edificios, acerad…  0   
##  2 1      1              dosel                                             89.9 
##  3 1      1              edificación erguida                                0   
##  4 1      1              suelo, herbáceas, no edificado ni cubierto        10.1 
##  5 2      1              construido, mobiliario (bordes edificios, acerad… 18.5 
##  6 2      1              dosel                                             80.2 
##  7 2      1              edificación erguida                                0   
##  8 2      1              suelo, herbáceas, no edificado ni cubierto         1.34
##  9 3      2              construido, mobiliario (bordes edificios, acerad…  0   
## 10 3      2              dosel                                             99.9 
## # … with 98 more rows

m_amb_ward_k2_tw <- m_amb_ward_k2 %>%
  group_by(variable) %>%
  summarise(
    p_valor_t = tryCatch(t.test(valor ~ grupos_ward_k2)$p.value, error = function(e) NA),
    p_valor_w = tryCatch(wilcox.test(valor ~ grupos_ward_k2)$p.value, error = function(e) NA)
    ) %>%
  arrange(p_valor_t)
m_amb_ward_k2_tw %>%
  kable(booktabs=T) %>%
  kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
  gsub(' NA ', '', .)

variable	p_valor_t	p_valor_w
construido, mobiliario (bordes edificios, acerado, bancos, postes…)	0.0405413	0.1910607
dosel	0.0696707	0.0743239
suelo, herbáceas, no edificado ni cubierto	0.6704152	0.8405846
edificación erguida	0.7283481	0.8108505

m_amb_ward_k2 %>% 
  group_by(variable) %>% 
        ggplot() + aes(x = grupos_ward_k2, y = valor, fill = grupos_ward_k2) + 
        geom_boxplot(lwd = 0.2) + 
        scale_fill_brewer(palette = 'Set1') +
        theme_bw(base_size=12) +
        theme(legend.position="none") +
        facet_wrap(~ variable, scales = 'free_y', ncol = 2,
                   labeller = label_wrap_gen(width = 35,multi_line = T))

# CON WARD K3 NO APARECEN EFECTOS SIGNIFICATIVOS. MEJOR USAR SOLO AGRUPAMIENTO WARD K2

En resumen, las variables de porcentaje de “construido…” y “dosel…” son las unicas que presentaron efectos significativos entre grupos, por lo que consideramos que estas son las que, con mayor probabilidad, se asocian con la composicion de especies diferenciada entre grupos. El grupo 1 se caracteriza por tener areas construidas predominantemente, y el 2 por tener predominio de dosel.

Realizar analisis de asociacion propiamente (correlacion punto biserial, funcion r.g)

phi_ward_k2 <- multipatt(
  mi_fam,
  grupos_ward_k2,
  func = "r.g",
  max.order = 1,
  control = how(nperm = 999))
summary(phi_ward_k2)

## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 11
##  Selected number of species: 3 
##  Number of species associated to 1 group: 3 
## 
##  List of species associated to each combination: 
## 
##  Group 1  #sps.  1 
##                        stat p.value    
## Dorymyrmex antillanus 0.749   0.001 ***
## 
##  Group 2  #sps.  2 
##                         stat p.value   
## Wasmannia auropunctata 0.655   0.002 **
## Cyphomyrmex rimosus    0.492   0.016 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ma_sf_grupos_ward_k2 <- ma_sf %>%
  inner_join(m_amb_ward_k2 %>%
               select(ENLACE=hex_id, grupos_ward_k2) %>%
               distinct() %>% 
               mutate(ENLACE=as.integer(ENLACE)))
st_geometry(ma_sf_grupos_ward_k2) <- 'geometry'
mapa_ma_sf_grupos_ward_k2 <- ma_sf_grupos_ward_k2 %>% 
  ggplot + aes(fill = grupos_ward_k2) +
  geom_sf()
if(interactive()) dev.new()
mapa_ma_sf_grupos_ward_k2

mapa_ma_sf <- ma_sf_grupos_ward_k2 %>% 
  select(-grupos_ward_k2) %>% 
  pivot_longer(-c(geometry, ENLACE)) %>% 
  group_by(name) %>% 
  mutate(value=scale(value)) %>% 
  ggplot + aes(fill = value) +
  scale_fill_viridis_b() +
  geom_sf() +
  facet_wrap(~name, ncol = 2)
if(interactive()) dev.new()
mapa_ma_sf

Análisis de ordenación

Analisis de ordenacion. No aporta patrones de interes, puesto que la inercia restringida es muy pequena (es habitual cuando se trabaja con datos de presencia/ausencia). Analisis descartado en este caso.

# mi_fam_t_rda <- rda(mi_fam_t ~ .,
#                     ma[-17,] %>% select(all_of(matches('construido|dosel'))),
#                     scale = T)
# summary(mi_fam_t_rda)
# RsquareAdj(mi_fam_t_rda)$adj.r.squared
# vif.cca(mi_fam_t_rda)
# escalado <- 1
# plot(mi_fam_t_rda,
#      scaling = escalado,
#      display = c("sp", "lc", "cn"),
#      main = paste("Triplot de RDA especies ~ var. GSL + ESA, escalamiento", escalado)
# )
# mi_fam_t_rda_sc1 <- scores(mi_fam_t_rda,
#          choices = 1:2,
#          scaling = escalado,
#          display = "sp"
#   )
# # text(mi_fam_t_rda, "species", col="red", cex=0.8, scaling=escalado)
# arrows(0, 0,
#        mi_fam_t_rda_sc1[, 1] * 0.9,
#        mi_fam_t_rda_sc1[, 2] * 0.9,
#        length = 0,
#        lty = 1,
#        col = "red"
# )

Diversidad, variabilidad estacional

Estimacion de riqueza (Chao)

Matriz combinada, ambas estaciones

Pool de especies, estimadores tradicionales de la riqueza (Chao, Jack1, Jack2, etc.). Esta matriz tiene muchas especies raras, por lo que los estimadores penalizan la riqueza capturada. Se considerara usar una matriz ponderada por el numero de nidos.

# mc_orig <- mc_todas$Combinada
specpool(mc_orig)

##     Species     chao  chao.se    jack1 jack1.se    jack2    boot  boot.se  n
## All      16 28.05357 16.54004 20.82143 2.564723 24.57011 18.0219 1.375846 28

specpool(mc_orig)[[1]]/specpool(mc_orig)[-c(3,5,8)]*100

##     Species     chao    jack1    jack2     boot        n
## All     100 57.03374 76.84391 65.11978 88.78088 57.14286

# Nidos
estimateR(colSums(mc_nidos))

##     S.obs   S.chao1  se.chao1     S.ACE    se.ACE 
## 16.000000 26.000000 10.258674 20.465909  2.201645

estimateR(colSums(mc_nidos))[[1]]/estimateR(colSums(mc_nidos))[-c(3,5,8)]*100

##     S.obs   S.chao1     S.ACE 
## 100.00000  61.53846  78.17879

Estimadores

Matriz original combinada: el error “Error in if (var_iChao2 > 0) { : missing value where TRUE/FALSE needed” impide ejecutar la estimacion para la matriz combinada.

# Incidencia
ChaoSpecies(data.frame(V1 = c(nrow(mc_orig), as.integer(colSums(mc_orig)))),
            datatype = 'incidence_freq', k=1, conf=0.95)

# Nmero de nidos
ChaoSpecies(data.frame(V1 = c(nrow(mc_nidos), as.integer(colSums(mc_nidos)))),
            datatype = 'abundance', k=1, conf=0.95)

## Warning: In this case, it can't estimate the variance of Homogeneous estimation

## 
## (1) BASIC DATA INFORMATION:
## 
##                                          Variable Value
##     Sample size                                 n   283
##     Number of observed species                  D    17
##     Coverage estimate for entire dataset        C 0.982
##     CV for entire dataset                      CV 1.538
##     Cut-off point                               k     1
## 
##                                                       Variable Value
##     Number of observed individuals for rare group       n_rare     5
##     Number of observed species for rare group           D_rare     5
##     Estimate of the sample coverage for rare group      C_rare     0
##     Estimate of CV for rare group in ACE               CV_rare     0
##     Estimate of CV1 for rare group in ACE-1           CV1_rare     0
##     Number of observed individuals for abundant group   n_abun   278
##     Number of observed species for abundant group       D_abun    12
## 
## NULL
## 
## 
## (2) SPECIES RICHNESS ESTIMATORS TABLE:
## 
##                               Estimate   s.e. 95%Lower 95%Upper
##     Homogeneous Model              Inf     NA       NA       NA
##     Homogeneous (MLE)           17.000  1.123   17.000   21.705
##     Chao1 (Chao, 1984)          26.965 10.235   18.888   69.583
##     Chao1-bc                    26.965 10.235   18.888   69.583
##     iChao1 (Chiu et al. 2014)   29.465 10.607   19.934   69.955
##     ACE (Chao & Lee, 1992)      26.965 10.235   18.888   69.583
##     ACE-1 (Chao & Lee, 1992)    26.965 10.235   18.888   69.583
##     1st order jackknife         21.982  3.154   18.597   32.544
##     2nd order jackknife         26.947  5.453   20.642   44.166
## 
## 
## (3) DESCRIPTION OF ESTIMATORS/MODELS:
## 
## Homogeneous Model: This model assumes that all species have the same incidence or detection probabilities. See Eq. (3.2) of Lee and Chao (1994) or Eq. (12a) in Chao and Chiu (2016b).
## 
## Chao2 (Chao, 1987): This approach uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Chao (1987) or Eq. (11a) in Chao and Chiu (2016b).
##      
## Chao2-bc: A bias-corrected form for the Chao2 estimator; see Chao (2005).
##   
## iChao2: An improved Chao2 estimator; see Chiu et al. (2014).
## 
## ICE (Incidence-based Coverage Estimator): A non-parametric estimator originally proposed by Lee and Chao (1994) in the context of capture-recapture data analysis. The observed species are separated as frequent and infrequent species groups; only data in the infrequent group are used to estimate the number of undetected species. The estimated CV for species in the infrequent group characterizes the degree of heterogeneity among species incidence probabilities. See Eq. (12b) of Chao and Chiu (2016b), which is an improved version of Eq. (3.18) in Lee and Chao (1994). This model is also called Model(h) in capture-recapture literature where h denotes "heterogeneity".
## 
## ICE-1: A modified ICE for highly-heterogeneous cases.
## 
## 1st order jackknife: It uses the frequency of uniques to estimate the number of undetected species; see Burnham and Overton (1978).
## 
## 2nd order jackknife: It uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Burnham and Overton (1978).
## 
## 95% Confidence interval: A log-transformation is used for all estimators so that the lower bound of the resulting interval is at least the number of observed species. See Chao (1987).

Usando matriz combinada, incluyendo singletons

Para fines didacticos, representamos la curva de acumulacion de especies de la matriz combinada, donde se observa una pronunciada pendiente y una escasa capacidad de alcanzarse la asintota. Decimos “para fines didacticos”, porque complementa el analisis posterior por grupos segun Ward.

mc_orig_lista <- sapply(rownames(mc_orig), function(x) as.numeric(mc_orig[x,]), simplify = F)
nasin_raref_pooled <- iNEXT::iNEXT(
  x = colSums(mc_orig), #if(is.list(mc_orig)) mc_orig_lista else mc_orig,
  q=0,
  knots = 400,
  datatype = 'incidence_freq')
acumulacion_especies_pooled <- iNEXT::ggiNEXT(nasin_raref_pooled, type=1) +
  theme_bw() +
  theme(
    text = element_text(size = 20),
    panel.background = element_rect(fill = 'white', colour = 'black'),
    panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
  ) +
  ylab('Riqueza de especies') +
  xlab('Número de sitios') +
  scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
  scale_color_manual(values = brewer.pal(8, 'Set2')) +
  scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies_pooled

Numero de nidos

mc_nidos_lista <- sapply(rownames(mc_nidos), function(x) as.numeric(mc_nidos[x,]), simplify = F)
nasin_raref_nidos <- iNEXT::iNEXT(
  x = colSums(mc_nidos), #if(is.list(mc_nidos)) mc_nidos_lista else mc_nidos,
  q=0,
  knots = 400,
  datatype = 'abundance')
acumulacion_especies_nidos <- iNEXT::ggiNEXT(nasin_raref_nidos, type=1) +
  theme_bw() +
  theme(
    text = element_text(size = 20),
    panel.background = element_rect(fill = 'white', colour = 'black'),
    panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
  ) +
  ylab('Riqueza de especies') +
  xlab('Número de sitios') +
  scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
  scale_color_manual(values = brewer.pal(8, 'Set2')) +
  scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies_nidos

Numero de nidos reescalado

mc_nidos_reescalado_lista <- sapply(rownames(mc_nidos_reescalado), function(x) as.numeric(mc_nidos_reescalado[x,]), simplify = F)
nasin_raref_nidos_reescalado <- iNEXT::iNEXT(
  x = colSums(mc_nidos_reescalado), #if(is.list(mc_nidos)) mc_nidos_lista else mc_nidos,
  q=0,
  knots = 400,
  datatype = 'abundance')
acumulacion_especies_nidos_reescalado <- iNEXT::ggiNEXT(nasin_raref_nidos_reescalado, type=1) +
  theme_bw() +
  theme(
    text = element_text(size = 20),
    panel.background = element_rect(fill = 'white', colour = 'black'),
    panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
  ) +
  ylab('Riqueza de especies') +
  xlab('Número de sitios') +
  scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
  scale_color_manual(values = brewer.pal(8, 'Set2')) +
  scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies_nidos_reescalado

Usando matriz sin especies singletons, usada en analisis de agrupamiento

Se refiere a la matriz filtrada donde se eliminaron especies poco representadas

specpool(mi_fam)

##     Species chao chao.se jack1 jack1.se    jack2     boot   boot.se  n
## All      11   11       0    11        0 10.10969 11.21473 0.4439909 27

ChaoSpecies(data.frame(V1 = c(nrow(mi_fam), as.numeric(colSums(mi_fam)))),
            datatype = 'incidence_freq', k=2, conf=0.95)

## Warning: In this case, it can't estimate the variance of Homogeneous estimation 
## 
## Warning: In this case, it can't estimate the variance of Model(h) estimation 
## 
## Warning: In this case, it can't estimate the variance of Model(h)-1 estimation 
## 
## Warning: In this case, it can't estimate the variance of 1st-order-jackknife estimation 
## 
## Warning: In this case, it can't estimate the variance of 2nd-order-jackknife estimation 
## 
## Warning: In this case, it can't estimate the variance of Model(h) estimation 
## 
## Warning: In this case, it can't estimate the variance of Model(h)-1 estimation

## 
## (1) BASIC DATA INFORMATION:
## 
##                                          Variable Value
##     Number of observed species                  D    11
##     Number of sampling units                    T    27
##     Total number of incidences                  U   104
##     Coverage estimate for entire dataset        C     1
##     CV for entire dataset                      CV 0.759
## 
##                                                       Variable Value
##     Cut-off point                                            k     2
##     Total number of incidences in infrequent group    U_infreq     2
##     Number of observed species for infrequent group   D_infreq     1
##     Estimated sample coverage for infrequent group    C_infreq     1
##     Estimated CV for infrequent group in ICE         CV_infreq 0.196
##     Estimated CV1 for infrequent group in ICE-1     CV1_infreq 0.196
##     Number of observed species for frequent group       D_freq    10
## 
##                            Q1 Q2
##     Incidence freq. counts  0  1
## 
## 
## (2) SPECIES RICHNESS ESTIMATORS TABLE:
## 
##                               Estimate  s.e. 95%Lower 95%Upper
##     Homogeneous Model            11.00 0.457       11   12.169
##     Chao2 (Chao, 1987)           11.00 0.457       11   12.169
##     Chao2-bc                     11.00 0.457       11   12.169
##     iChao2 (Chiu et al. 2014)    11.00 0.457       11   12.169
##     ICE (Lee & Chao, 1994)       11.00 0.457       11   12.169
##     ICE-1 (Lee & Chao, 1994)     11.00 0.457       11   12.169
##     1st order jackknife          11.00 0.457       11   12.169
##     2nd order jackknife          10.11    NA       NA       NA
## 
## 
## (3) DESCRIPTION OF ESTIMATORS/MODELS:
## 
## Homogeneous Model: This model assumes that all species have the same incidence or detection probabilities. See Eq. (3.2) of Lee and Chao (1994) or Eq. (12a) in Chao and Chiu (2016b).
## 
## Chao2 (Chao, 1987): This approach uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Chao (1987) or Eq. (11a) in Chao and Chiu (2016b).
##      
## Chao2-bc: A bias-corrected form for the Chao2 estimator; see Chao (2005).
##   
## iChao2: An improved Chao2 estimator; see Chiu et al. (2014).
## 
## ICE (Incidence-based Coverage Estimator): A non-parametric estimator originally proposed by Lee and Chao (1994) in the context of capture-recapture data analysis. The observed species are separated as frequent and infrequent species groups; only data in the infrequent group are used to estimate the number of undetected species. The estimated CV for species in the infrequent group characterizes the degree of heterogeneity among species incidence probabilities. See Eq. (12b) of Chao and Chiu (2016b), which is an improved version of Eq. (3.18) in Lee and Chao (1994). This model is also called Model(h) in capture-recapture literature where h denotes "heterogeneity".
## 
## ICE-1: A modified ICE for highly-heterogeneous cases.
## 
## 1st order jackknife: It uses the frequency of uniques to estimate the number of undetected species; see Burnham and Overton (1978).
## 
## 2nd order jackknife: It uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Burnham and Overton (1978).
## 
## 95% Confidence interval: A log-transformation is used for all estimators so that the lower bound of the resulting interval is at least the number of observed species. See Chao (1987).

Curva de acumulacion de especies. Primero generamos la matriz combinada (pooled) por grupos del agrupamiento Ward K2.

mi_fam_ward <- mi_fam %>%
  mutate(g = grupos_ward_k2) %>%
  group_by(g) %>%
  summarise_all(sum) %>%
  select(-g) %>% 
  mutate(N = nrow(mi_fam)) %>% 
  relocate(N, .before = 1) %>% 
  data.frame
mi_fam_ward

##    N Pheidole.subarmata Dorymyrmex.antillanus Paratrechina.longicornis
## 1 27                 16                    16                        3
## 2 27                  9                     2                        5
##   Solenopsis.geminata Pheidole.jelskii Mycetomoellerius.jamaicensis
## 1                  11                6                            1
## 2                   9                2                            1
##   Tetramorium.lanuginosum Odontomachus.ruginodis Wasmannia.auropunctata
## 1                       1                      2                      0
## 2                       2                      3                      6
##   Cyphomyrmex.rimosus Pogonomyrmex.schmitti
## 1                   1                     1
## 2                   5                     2

Posteriormente, generamos la curva de acumlacion.

nasin_raref <- iNEXT::iNEXT(
  x = t(mi_fam_ward),
  q=0,
  knots = 400,
  datatype = 'incidence_freq')
acumulacion_especies <- iNEXT::ggiNEXT(nasin_raref, type=1) +
  theme_bw() +
  theme(
    text = element_text(size = 20),
    panel.background = element_rect(fill = 'white', colour = 'black'),
    panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
  ) +
  ylab('Riqueza de especies') +
  xlab('Número de sitios') +
  scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
  scale_color_manual(values = brewer.pal(8, 'Set2')) +
  scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies

Por estaciones

Temporada seca

#' ### Número de sitios, tanto en matriz de comunidad como en ambiental
#' Verifica que coinciden
nrow(mc_todas$Marzo) #En la matriz de comunidad

## [1] 27

nrow(ma) #En la matriz ambiental

## [1] 28

#' ### Riqueza numérica de especies (usando matriz de comunidad) por quadrat
#' Nota: cargar paquete vegan arriba, en el área de paquetes
specnumber(mc_todas$Marzo)

##  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 
##  2  1  3  3  2  2  2  2  1  2  5  5  1  2  2  1  3  3  3  4  2  1  2  4  2  2 
## 28 
##  1

sort(specnumber(mc_todas$Marzo)) # Ordenados ascendentemente

##  3 10 14 17 23 28  2  6  7  8  9 11 15 16 22 24 26 27  4  5 18 19 20 21 25 12 
##  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  3  3  3  3  3  4  4  5 
## 13 
##  5

summary(specnumber(mc_todas$Marzo)) # Resumen estadístico

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.000   2.000   2.333   3.000   5.000

#' ### Riqueza numérica de toda la "comunidad"
specnumber(colSums(mc_todas$Marzo))

## [1] 12

Comparar matrices entre estaciones con TBI

# TBI
mc_marzo_para_tbi <- cbind(
  rbind(rownames(mc_pooled)[!rownames(mc_pooled) %in% rownames(mc_marzo)], mc_marzo) %>% sapply(., as.integer),
  mc_mayo[,!colnames(mc_mayo) %in% colnames(mc_marzo)] %>% mutate_all(.funs = ~ifelse(.>0, 0, 0)))
colnames(mc_marzo)[!colnames(mc_marzo) %in% colnames(mc_mayo)] #Presentes en marzo, ausentes en mayo

## [1] "Pheidole moerens"        "Tetramorium bicarinatum"

colnames(mc_mayo)[!colnames(mc_mayo) %in% colnames(mc_marzo)] #Presentes en mayo, ausentes en marzo

## [1] "Brachymyrmex heeri"    "Pogonomyrmex schmitti" "Monomorium floricola" 
## [4] "Acropyga dubitata"

mc_mayo_para_tbi <- cbind(
  mc_mayo,
  mc_marzo_para_tbi[,!colnames(mc_marzo_para_tbi) %in% colnames(mc_mayo)] %>% mutate_all(.funs = ~ifelse(.>0, 0, 0)))
mc_mayo_para_tbi <- mc_mayo_para_tbi[,colnames(mc_marzo_para_tbi)]
all.equal(colnames(mc_marzo_para_tbi), colnames(mc_mayo_para_tbi))

## [1] TRUE

rownames(mc_marzo_para_tbi) == rownames(mc_mayo_para_tbi)

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

tbi_porc_diferencias <- TBI(mc_marzo_para_tbi, mc_mayo_para_tbi, method = '%difference', nperm = 999, test.t.perm = T)
tbi_porc_diferencias

## $TBI
##  [1] 0.7142857 0.5000000 1.0000000 0.3333333 0.3333333 0.5555556 1.0000000
##  [8] 0.5000000 0.6000000 0.3333333 0.5000000 0.4285714 0.3333333 1.0000000
## [15] 0.5000000 0.5000000 0.0000000 0.2000000 0.6000000 0.3333333 0.3333333
## [22] 0.2000000 1.0000000 0.6000000 0.4285714 0.0000000 0.3333333 0.6000000
## 
## $p.TBI
##  [1] 0.244 0.592 0.153 0.820 0.803 0.428 0.175 0.576 0.404 0.828 0.574 0.671
## [13] 0.820 0.171 0.568 0.538 0.997 0.912 0.412 0.828 0.823 0.931 0.173 0.382
## [25] 0.645 1.000 0.805 0.397
## 
## $p.adj
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $BCD.mat
##         B/(2A+B+C) C/(2A+B+C) D=(B+C)/(2A+B+C) Change
## Site.1   0.7142857  0.0000000        0.7142857    -  
## Site.2   0.2500000  0.2500000        0.5000000    0  
## Site.3   0.5000000  0.5000000        1.0000000    0  
## Site.4   0.1666667  0.1666667        0.3333333    0  
## Site.5   0.1666667  0.1666667        0.3333333    0  
## Site.6   0.0000000  0.5555556        0.5555556    +  
## Site.7   0.2222222  0.7777778        1.0000000    +  
## Site.8   0.2500000  0.2500000        0.5000000    0  
## Site.9   0.2000000  0.4000000        0.6000000    +  
## Site.10  0.0000000  0.3333333        0.3333333    +  
## Site.11  0.0000000  0.5000000        0.5000000    +  
## Site.12  0.4285714  0.0000000        0.4285714    -  
## Site.13  0.2222222  0.1111111        0.3333333    -  
## Site.14  0.5000000  0.5000000        1.0000000    0  
## Site.15  0.0000000  0.5000000        0.5000000    +  
## Site.16  0.2500000  0.2500000        0.5000000    0  
## Site.17  0.0000000  0.0000000        0.0000000    0  
## Site.18  0.2000000  0.0000000        0.2000000    -  
## Site.19  0.4000000  0.2000000        0.6000000    -  
## Site.20  0.1666667  0.1666667        0.3333333    0  
## Site.21  0.1111111  0.2222222        0.3333333    +  
## Site.22  0.0000000  0.2000000        0.2000000    +  
## Site.23  0.5000000  0.5000000        1.0000000    0  
## Site.24  0.2000000  0.4000000        0.6000000    +  
## Site.25  0.2857143  0.1428571        0.4285714    -  
## Site.26  0.0000000  0.0000000        0.0000000    0  
## Site.27  0.0000000  0.3333333        0.3333333    +  
## Site.28  0.0000000  0.6000000        0.6000000    +  
## 
## $BCD.summary
##  mean(B/den) mean(C/den)   mean(D)   B/(B+C)   C/(B+C) Change
##    0.2047902   0.2866497 0.4914399 0.4167147 0.5832853    +  
## 
## $t.test_B.C
##                 mean(C-B)     Stat   p.param p.perm   p<=0.05
## Paired t.test  0.08185941 1.437708 0.1620069  0.174          
## 
## $BC
## [1] NA
## 
## attr(,"class")
## [1] "TBI"

tbi_jaccard <- TBI(mc_marzo_para_tbi, mc_mayo_para_tbi, method = 'jaccard', nperm = 999, test.t.perm = T)
tbi_jaccard

## $TBI
##  [1] 0.8333333 0.6666667 1.0000000 0.5000000 0.5000000 0.7142857 1.0000000
##  [8] 0.6666667 0.7500000 0.5000000 0.6666667 0.6000000 0.5000000 1.0000000
## [15] 0.6666667 0.6666667 0.0000000 0.3333333 0.7500000 0.5000000 0.5000000
## [22] 0.3333333 1.0000000 0.7500000 0.6000000 0.0000000 0.5000000 0.7500000
## 
## $p.TBI
##  [1] 0.210 0.582 0.171 0.818 0.806 0.395 0.166 0.553 0.397 0.820 0.580 0.656
## [13] 0.795 0.170 0.556 0.554 0.999 0.919 0.435 0.817 0.829 0.923 0.163 0.398
## [25] 0.644 0.999 0.825 0.400
## 
## $p.adj
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $BCD.mat
##         B/(A+B+C) C/(A+B+C) D=(B+C)/(A+B+C) Change
## Site.1  0.8333333 0.0000000       0.8333333    -  
## Site.2  0.3333333 0.3333333       0.6666667    0  
## Site.3  0.5000000 0.5000000       1.0000000    0  
## Site.4  0.2500000 0.2500000       0.5000000    0  
## Site.5  0.2500000 0.2500000       0.5000000    0  
## Site.6  0.0000000 0.7142857       0.7142857    +  
## Site.7  0.2222222 0.7777778       1.0000000    +  
## Site.8  0.3333333 0.3333333       0.6666667    0  
## Site.9  0.2500000 0.5000000       0.7500000    +  
## Site.10 0.0000000 0.5000000       0.5000000    +  
## Site.11 0.0000000 0.6666667       0.6666667    +  
## Site.12 0.6000000 0.0000000       0.6000000    -  
## Site.13 0.3333333 0.1666667       0.5000000    -  
## Site.14 0.5000000 0.5000000       1.0000000    0  
## Site.15 0.0000000 0.6666667       0.6666667    +  
## Site.16 0.3333333 0.3333333       0.6666667    0  
## Site.17 0.0000000 0.0000000       0.0000000    0  
## Site.18 0.3333333 0.0000000       0.3333333    -  
## Site.19 0.5000000 0.2500000       0.7500000    -  
## Site.20 0.2500000 0.2500000       0.5000000    0  
## Site.21 0.1666667 0.3333333       0.5000000    +  
## Site.22 0.0000000 0.3333333       0.3333333    +  
## Site.23 0.5000000 0.5000000       1.0000000    0  
## Site.24 0.2500000 0.5000000       0.7500000    +  
## Site.25 0.4000000 0.2000000       0.6000000    -  
## Site.26 0.0000000 0.0000000       0.0000000    0  
## Site.27 0.0000000 0.5000000       0.5000000    +  
## Site.28 0.0000000 0.7500000       0.7500000    +  
## 
## $BCD.summary
##  mean(B/den) mean(C/den)   mean(D)   B/(B+C)   C/(B+C) Change
##    0.2549603   0.3610261 0.6159864 0.4139058 0.5860942    +  
## 
## $t.test_B.C
##                 mean(C-B)     Stat   p.param p.perm   p<=0.05
## Paired t.test   0.1060658 1.455444 0.1570758  0.152          
## 
## $BC
## [1] NA
## 
## attr(,"class")
## [1] "TBI"

tbi_sorensen <- TBI(mc_marzo_para_tbi, mc_mayo_para_tbi, method = 'sorensen', nperm = 999, test.t.perm = T)
tbi_sorensen

## $TBI
##  [1] 0.7142857 0.5000000 1.0000000 0.3333333 0.3333333 0.5555556 1.0000000
##  [8] 0.5000000 0.6000000 0.3333333 0.5000000 0.4285714 0.3333333 1.0000000
## [15] 0.5000000 0.5000000 0.0000000 0.2000000 0.6000000 0.3333333 0.3333333
## [22] 0.2000000 1.0000000 0.6000000 0.4285714 0.0000000 0.3333333 0.6000000
## 
## $p.TBI
##  [1] 0.219 0.580 0.170 0.805 0.814 0.420 0.157 0.564 0.415 0.819 0.554 0.643
## [13] 0.814 0.167 0.560 0.565 1.000 0.940 0.435 0.827 0.830 0.923 0.148 0.410
## [25] 0.622 1.000 0.789 0.400
## 
## $p.adj
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $BCD.mat
##         B/(2A+B+C) C/(2A+B+C) D=(B+C)/(2A+B+C) Change
## Site.1   0.7142857  0.0000000        0.7142857    -  
## Site.2   0.2500000  0.2500000        0.5000000    0  
## Site.3   0.5000000  0.5000000        1.0000000    0  
## Site.4   0.1666667  0.1666667        0.3333333    0  
## Site.5   0.1666667  0.1666667        0.3333333    0  
## Site.6   0.0000000  0.5555556        0.5555556    +  
## Site.7   0.2222222  0.7777778        1.0000000    +  
## Site.8   0.2500000  0.2500000        0.5000000    0  
## Site.9   0.2000000  0.4000000        0.6000000    +  
## Site.10  0.0000000  0.3333333        0.3333333    +  
## Site.11  0.0000000  0.5000000        0.5000000    +  
## Site.12  0.4285714  0.0000000        0.4285714    -  
## Site.13  0.2222222  0.1111111        0.3333333    -  
## Site.14  0.5000000  0.5000000        1.0000000    0  
## Site.15  0.0000000  0.5000000        0.5000000    +  
## Site.16  0.2500000  0.2500000        0.5000000    0  
## Site.17  0.0000000  0.0000000        0.0000000    0  
## Site.18  0.2000000  0.0000000        0.2000000    -  
## Site.19  0.4000000  0.2000000        0.6000000    -  
## Site.20  0.1666667  0.1666667        0.3333333    0  
## Site.21  0.1111111  0.2222222        0.3333333    +  
## Site.22  0.0000000  0.2000000        0.2000000    +  
## Site.23  0.5000000  0.5000000        1.0000000    0  
## Site.24  0.2000000  0.4000000        0.6000000    +  
## Site.25  0.2857143  0.1428571        0.4285714    -  
## Site.26  0.0000000  0.0000000        0.0000000    0  
## Site.27  0.0000000  0.3333333        0.3333333    +  
## Site.28  0.0000000  0.6000000        0.6000000    +  
## 
## $BCD.summary
##  mean(B/den) mean(C/den)   mean(D)   B/(B+C)   C/(B+C) Change
##    0.2047902   0.2866497 0.4914399 0.4167147 0.5832853    +  
## 
## $t.test_B.C
##                 mean(C-B)     Stat   p.param p.perm   p<=0.05
## Paired t.test  0.08185941 1.437708 0.1620069  0.151          
## 
## $BC
## [1] NA
## 
## attr(,"class")
## [1] "TBI"

# Calculando Jaccard de forma independiente
betadiver(rbind(colSums(mc_marzo_para_tbi), colSums(mc_mayo_para_tbi)), method = 'j')

##       1
## 2 0.625

# Para contrastar con tabla de valores significativos.

Resumen de resultado sobre TBI y diversidad Beta estacional

Calculando distintos indices de disimilaridad (%diferencias, Jaccard, Sorensen), y aplicando pruebas de permutacion a la diversidad Beta temporal (T1 y T2), se obtiene un resultado consistente sobre un cambio positivo en la diversidad de especies (ganancia de especies), pero este cambio no es significativo en terminos estadisticos.

Asimismo, tras calcular el indice de Jaccard mediante betadiver, contrastamos con la tabla de valores significativos de Raymundo Real, evidenciamos que existen diferencias significativas en la diversidad Beta comparando las matrices combinadas (pooled) de las estaciones para un nivel de significancia de 0.05, pero no para un nivel de significancia de 0.01. Esto sugiere que las diferencias en cuanto a la composicion de hormigas nidificantes entre ambas estaciones son realmente bajas.

Analisis univariado de la precipitacion, meses marzo y mayo de 2022

precip_orig <- read_excel('JARDIN BOTANICO (JOSE RAMON).xls', range = "A9:I40")
precip <- precip_orig %>% 
  mutate_all(~ gsub('INAP', '0', .x)) %>%
  mutate_all(~ gsub('-', NA, .x)) %>% 
  mutate_all(~ as.numeric(.x))
sapply(precip, summary)

## $DIA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0     8.5    16.0    16.0    23.5    31.0 
## 
## $ENE
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2484  0.0000  3.3000 
## 
## $FEB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   1.518   1.925  11.400       3 
## 
## $MAR
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   7.116  13.200  34.800 
## 
## $ABR
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   7.323   6.025  63.200       1 
## 
## $MAY
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   3.035   2.250  29.500 
## 
## $JUN
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   3.843   1.300  68.800       1 
## 
## $JUL
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   9.097  13.950  47.800 
## 
## $AGO
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   7.353  17.325  34.500       1

sapply(precip, sum, na.rm=T)

##   DIA   ENE   FEB   MAR   ABR   MAY   JUN   JUL   AGO 
## 496.0   7.7  42.5 220.6 219.7  94.1 115.3 282.0 220.6

precip_mar_colec <- precip %>% 
  select(DIA, MAR) %>% 
  filter(DIA %in% c(21:25, 30, 31)) %>% 
  pull(MAR)
sum(precip_mar_colec)

## [1] 27.1

precip_may_colec <- precip %>% 
  select(DIA, MAY) %>% 
  filter(DIA %in% c(10:13, 16:19, 31)) %>% 
  pull(MAY)
sum(precip_may_colec)

## [1] 46.7

Ecologa espacial

Generar hexagonos topologicamente correctos

Autocorrelación espacial mediante prueba Mantel (matrices de distancia)

library(sp)
# Matriz transformada Hellinger
mc_nidos_hell <- decostand(mc_nidos, 'hellinger')
# Matriz ambiental con centroides
st_geometry(hex_inters) <- 'geom'
ma_sf <- hex_inters %>% 
  inner_join(ma %>%
               rownames_to_column('ENLACE') %>%
               mutate(ENLACE = as.integer(ENLACE)))
ma_sp <- ma_sf %>% as_Spatial()
centroides <- ma_sf %>% st_centroid
ma_xy <- centroides %>% st_coordinates %>% as.data.frame
(vecindad <- centroides %>% dnearneigh(d1 = 0, d2 = 75))

## Neighbour list object:
## Number of regions: 28 
## Number of nonzero links: 128 
## Percentage nonzero weights: 16.32653 
## Average number of links: 4.571429

plot(ma_sp)
plot(vecindad, coords = ma_xy, add=T, col = 'red')

(pesos_b <- nb2listw(vecindad, style = 'B'))

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 28 
## Number of nonzero links: 128 
## Percentage nonzero weights: 16.32653 
## Average number of links: 4.571429 
## 
## Weights style: B 
## Weights constants summary:
##    n  nn  S0  S1   S2
## B 28 784 128 256 2544

# Matriz sin tendencia
mc_nidos_sin_tendencia <- resid(
  lm(as.matrix(mc_nidos_hell) ~ .,
     data = ma_xy))
(autocor_global_residuos_mc <- sapply(
  dimnames(mc_nidos_sin_tendencia)[[2]],
  function(x)
    moran.mc(
      x = mc_nidos_sin_tendencia[,x],
      listw = pesos_b,
      zero.policy = T,
      nsim = 9999),
    simplify = F))

## $`Pheidole subarmata`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.12312, observed rank = 2306, p-value = 0.7694
## alternative hypothesis: greater
## 
## 
## $`Dorymyrmex antillanus`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.084899, observed rank = 3532, p-value = 0.6468
## alternative hypothesis: greater
## 
## 
## $`Paratrechina longicornis`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.13807, observed rank = 1794, p-value = 0.8206
## alternative hypothesis: greater
## 
## 
## $`Solenopsis geminata`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = 0.14283, observed rank = 9320, p-value = 0.068
## alternative hypothesis: greater
## 
## 
## $`Pheidole jelskii`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.16881, observed rank = 1021, p-value = 0.8979
## alternative hypothesis: greater
## 
## 
## $`Pheidole moerens`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.093401, observed rank = 1420, p-value = 0.858
## alternative hypothesis: greater
## 
## 
## $`Mycetomoellerius jamaicensis`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.082011, observed rank = 3406, p-value = 0.6594
## alternative hypothesis: greater
## 
## 
## $`Tetramorium lanuginosum`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.093645, observed rank = 2933, p-value = 0.7067
## alternative hypothesis: greater
## 
## 
## $`Odontomachus ruginodis`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.29487, observed rank = 30, p-value = 0.997
## alternative hypothesis: greater
## 
## 
## $`Wasmannia auropunctata`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.093884, observed rank = 3161, p-value = 0.6839
## alternative hypothesis: greater
## 
## 
## $`Cyphomyrmex rimosus`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = 0.024807, observed rank = 7402, p-value = 0.2598
## alternative hypothesis: greater
## 
## 
## $`Tetramorium bicarinatum`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.097105, observed rank = 1337, p-value = 0.8663
## alternative hypothesis: greater
## 
## 
## $`Brachymyrmex heeri`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.096045, observed rank = 287, p-value = 0.9713
## alternative hypothesis: greater
## 
## 
## $`Pogonomyrmex schmitti`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.19293, observed rank = 531, p-value = 0.9469
## alternative hypothesis: greater
## 
## 
## $`Monomorium floricola`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = 0.018724, observed rank = 8068, p-value = 0.1932
## alternative hypothesis: greater
## 
## 
## $`Acropyga dubitata`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = 0.018724, observed rank = 8091, p-value = 0.1909
## alternative hypothesis: greater

(autocor_global_residuos <- sapply(
  dimnames(mc_nidos_sin_tendencia)[[2]],
  function(x)
    moran.test(
      x = mc_nidos_sin_tendencia[,x],
      listw = pesos_b,
      zero.policy = T),
    simplify = F))

## $`Pheidole subarmata`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.77838, p-value = 0.7818
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.12312262       -0.03703704        0.01223152 
## 
## 
## $`Dorymyrmex antillanus`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.41818, p-value = 0.6621
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.08489939       -0.03703704        0.01309983 
## 
## 
## $`Paratrechina longicornis`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.92421, p-value = 0.8223
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.13807459       -0.03703704        0.01195157 
## 
## 
## $`Solenopsis geminata`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = 1.576, p-value = 0.05752
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.14282977       -0.03703704        0.01302614 
## 
## 
## $`Pheidole jelskii`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -1.1842, p-value = 0.8818
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.16880693       -0.03703704        0.01238092 
## 
## 
## $`Pheidole moerens`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -1.0731, p-value = 0.8584
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.093401039      -0.037037037       0.002759073 
## 
## 
## $`Mycetomoellerius jamaicensis`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.50945, p-value = 0.6948
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.082011111      -0.037037037       0.007793372 
## 
## 
## $`Tetramorium lanuginosum`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.62151, p-value = 0.7329
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.093645079      -0.037037037       0.008295785 
## 
## 
## $`Odontomachus ruginodis`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -2.3369, p-value = 0.9903
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.29487015       -0.03703704        0.01217280 
## 
## 
## $`Wasmannia auropunctata`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.50626, p-value = 0.6937
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.09388384       -0.03703704        0.01260874 
## 
## 
## $`Cyphomyrmex rimosus`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = 0.56388, p-value = 0.2864
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.02480714       -0.03703704        0.01202887 
## 
## 
## $`Tetramorium bicarinatum`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -1.1375, p-value = 0.8723
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.097105008      -0.037037037       0.002788338 
## 
## 
## $`Brachymyrmex heeri`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -1.904, p-value = 0.9715
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     -0.0960451331     -0.0370370370      0.0009604943 
## 
## 
## $`Pogonomyrmex schmitti`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -1.4674, p-value = 0.9289
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.19292906       -0.03703704        0.01128672 
## 
## 
## $`Monomorium floricola`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = 0.88658, p-value = 0.1877
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.018724330      -0.037037037       0.003955756 
## 
## 
## $`Acropyga dubitata`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = 0.88658, p-value = 0.1877
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.018724330      -0.037037037       0.003955756

pvalor_moran <- sapply(names(autocor_global_residuos),
                       function(x) autocor_global_residuos[[x]]$p.value)
tabla_patrones_espaciales <- data.frame(
  Especie = names(pvalor_moran),
  ValorP = pvalor_moran,
  Patron = ifelse(pvalor_moran<0.1, 'Agregado', 'Aleatorio'))
tabla_patrones_espaciales %>% kable

	Especie	ValorP	Patron
Pheidole subarmata	Pheidole subarmata	0.7818267	Aleatorio
Dorymyrmex antillanus	Dorymyrmex antillanus	0.6620916	Aleatorio
Paratrechina longicornis	Paratrechina longicornis	0.8223113	Aleatorio
Solenopsis geminata	Solenopsis geminata	0.0575183	Agregado
Pheidole jelskii	Pheidole jelskii	0.8818410	Aleatorio
Pheidole moerens	Pheidole moerens	0.8583758	Aleatorio
Mycetomoellerius jamaicensis	Mycetomoellerius jamaicensis	0.6947807	Aleatorio
Tetramorium lanuginosum	Tetramorium lanuginosum	0.7328685	Aleatorio
Odontomachus ruginodis	Odontomachus ruginodis	0.9902783	Aleatorio
Wasmannia auropunctata	Wasmannia auropunctata	0.6936615	Aleatorio
Cyphomyrmex rimosus	Cyphomyrmex rimosus	0.2864181	Aleatorio
Tetramorium bicarinatum	Tetramorium bicarinatum	0.8723456	Aleatorio
Brachymyrmex heeri	Brachymyrmex heeri	0.9715441	Aleatorio
Pogonomyrmex schmitti	Pogonomyrmex schmitti	0.9288624	Aleatorio
Monomorium floricola	Monomorium floricola	0.1876521	Aleatorio
Acropyga dubitata	Acropyga dubitata	0.1876521	Aleatorio

# mc_nidos_sin_tendencia_d <- dist(mc_nidos_sin_tendencia)
# (mc_nidos_correlograma <- mantel.correlog(
#   mc_nidos_sin_tendencia_d,
#   XY = ma_xy,
#   nperm = 999))
# plot(mc_nidos_correlograma)